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Abstract. A recent third-order, essentially non-oscillatory 
central scheme to advance the equations of single-fluid mag- 
netohydrodynamics (MHD) in time has been implemented 
into a new numerical code. This code operates on a 3-D 
Cartesian, non-staggered grid, and is able to handle shock- 
like gradients without producing spurious oscillations. 
To demonstrate the suitability of our code for the simulation 
of coronal mass ejections (CMEs) and similar heliospheric 
transients, we present selected results from test cases and 
perform studies of the solar wind expansion during phases of 
minimum solar activity. We can demonstrate convergence of 
the system into a stable Parker-like steady state for both hy- 
drodynamic and MHD winds. The model is subsequently ap- 
plied to expansion studies of CME-like plasma bubbles, and 
their evolution is monitored until a stationary state similar to 
the initial one is achieved. In spite of the model's (current) 
simplicity, we can confirm the CME's nearly self- similar 
evolution close to the Sun, thus highlighting the importance 
of detailed modelling especially at small heliospheric radii. 
Additionally, alternative methods to implement boundary 
conditions at the coronal base, as well as strategies to ensure 
a solenoidal magnetic field, are discussed and evaluated. 

Keywords. Interplanetary physics (Solar wind plasma) - 
Solar physics, astrophysics, and astronomy (Flares and mass 
ejections) - Space plasma physics (Numerical simulation 
studies) 



1 Introduction 

Coronal mass ejections (CMEs) moved into the focus of 
several research activities during recent years. Besides a 
variety of observational data resulting from SOHO ( jPickj 



elaT) [2006]), SMEI (Webb et al.| 12006] ), and, very recently, 
STEREO ( Vourlidas et al. ~ |2Q07| ), significant progress has 
also been achieved with the numerical modelling of CMEs, 
see, e.g., the reviews by Aschwanden et al. ( 2006]) and |Forbes 



|et al.| ( [2006] ). The motivation for the various activities is at 
least fourfold. First, CMEs are amongst the main mediators 
of the influence of the Sun on the inner heliosphere, partic- 
ularly on the Earth and its environment, where they signifi- 
cantly co-determine the space weather conditions. There is, 
in view of ever advancing technology that is increasingly sen- 
sitive — if not vulnerable — to space weather effects, strong 
interest in an understanding of the latter. An even stronger 
driver of research activities is provided with the recognition 
that space weather phenomena offer valuable opportunities 
to study many aspects of plasma astrophysics in great detail 



( Sch erer et aL] [2005] |Bothmer a nd Daglfs] [2006] |Schwenn| 
2006|. Second, as a consequence of the shocks driven by 
CMEs, they serve as particle accelerators that do not only 
contribute to space weather effects, but can be used to study 



1999; Mewaldt 



the actual acceleration processes ( [Reames 
et al.| |2005| |Li et al.| |2005| ), which are expected to oc- 
cur in other astrophysical systems as well (Eichler, 2006). 
Third, with the recent launch of the two- spacecraft mission 
STEREO ( [Kaiser! [2005] ), the full three-dimensional struc- 
ture of CMEs can be observed both remotely and in- situ for 
the first time. First results have already been reported by, 
e.g., |Howard and Tappin] ( [2008] ) and |Vourlidas et al.| ( |2007] ). 
And, fourth, the magnetohydrodynamic (MHD) modelling 
of CMEs provides an excellent testbed for numerical codes. 
Although not strongly motivated by CME physics, this is ac- 
tually one of the main drivers of model development, as is 
manifest with numerous approaches documented in the lit- 
erature. These various approaches can be ordered into three 
groups. There is (i) principal modelling that is either analyt- 
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extended corona, i.e. a few tens of solar radii (e.g., Mikic 
and Linker] 1994), and (iii) global modelling covering the 
inner heliosphere from the solar surface out to 1 AU and be- 
yond (e.g.,|Manchester et al.[|2004[|Odstrcil et al. [12005} |T6th 



where 



[eTaIl|2005HRiley et al.[|2006] r 
Despite these intensified efforts and activity regarding the 
study of CMEs, there remains both a number of unsolved 
problems and various modelling deficiencies. For example, 
the acceleration and heating processes of the plasma near the 
coronal base are — even nearly 50 years after the 'discovery' 
of the solar wind — still not known (Cranm er et al.| |2007| ), 
and there is also no agreement on the processes that actu- 
ally initiate CMEs ( [Forbes etliL) |2006| ). Also, their prop- 
agation and evolution in size and shape is by far not fully 
understood in all detail, and neither is their interaction with 
the background solar wind (Jac obs et al.[ |2007] ), with other 
CMEs ( Gopalswamy et al. 2001 ), and with planetary mag 
netospheres (e.g., |Groth et al.| |2000t |Ip and Koppl |2002 



Regarding the model formulations underlying the numerical 
simulation of CMEs, particularly the (non-thermal) heating 
of the plasma is mostly treated in a rather simplified manner 
via ad-hoc heating functions (e.g., |Groth et al.[ 2000; Manch- 
ester et al. 2004), variable adiabatic indices 7 = 7(1*) (e.g., 



Lugaz et al. 2007|), o r phenomenological heating functions 



(e.g., |Usmano v et al.[ 
|etal.| ( |2008| ). 



2000), for a discussion see Fichtner 



With the intention to address several of the above-mentioned 
problems, we have applied our recently developed CWENO- 
based MHD code ( [Kleimann et al.[ 2004|), which primor- 
dially originated from that by Grauer an d Marliani ( 2000 ), 
to the CME expansion problem. To our knowledge, this 
is the first published paper describing the application of a 
CWENO-based numerical code to an MHD problem related 
to space physics. 

In the following, we describe the model (Sect. [2]) and its nu- 
merical realization (Sect. [3] to [5]), present results of analyses 
of both the propagation of individual and the interaction of 
two CMEs, and suggest a possible connection of our find- 
ings to observations (Sect. [6]). 

2 Governing equations 

Choosing the normalization constants summarized in Tab. [T] 
the set of MHD equations for mass density p, flow velocity 
u, magnetic field strength B, and gas pressure p reads (in 
dimensionless form): 



d t e 



^ + V-(pn)=0 
d t (p u) + V - [puu 
+(p+\\B\\ 2 /2)i-BB] =pg 

d t B + V • (uB - Bu) = 

-V- [(e+p+||S|| 2 /2)n 

-(u-B)B]=p(Q- 



ug) 



(1) 

(2) 
(3) 

(4) 



r/r 2 e r 
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p \\u 



and 

IBII 2 



fp/(7-l)=7^1 
\ : 7 = 1 



respectively denote gravity (with 

T := (GM Q )/(R Q c 2 s ) = 11.49 



(5) 
(6) 

(7) 



in normalized units, cf. Tab. [I]), and the total energy density 
of a plasma with adiabatic exponent 7. Throughout this 
paper, || • || is used to denote the norm of a vector (i.e. 
||X|J_= \JX • X for any vector X), and the symbol I 
Eq. d2b denotes the unit tensor (i.e. (1)^- = Sij). 



m 



Table 1. Normalization constants used for Eqs. ([T}j4]). M©, Rq, 
/cb , m p , po , and c s denote the solar mass, the solar radius, the Boltz- 
mann constant, the proton rest mass, the electrical permeability of 
free space, and the isothermal sound speed, respectively. 



Quantity 


Normalization 


(solar) mass 


M 




2.0 x 10 30 kg 


length 


Lo := Rq 




7.0 x 10 8 m 


temperature T 


To 




1.0 x 10 6 K 


number density n = 


p/m p n 




1.0 x 10 14 m- 3 


plasma pressure p 


2 n k B To 




2.8 x 10" 3 Pa 


velocity u c s 


:= y/2 k B T /m p 




1.3 x 10 5 ms" 1 


time t 


Lo/cs 
m p n (c s ) 2 




5.5 x 10 3 s 


energy density e 




2.8 x 10" 3 Jm" 


mag. induction B 


c s yjjio m p no 




4.2 x 10" 5 T 


heating rate Q 


(c s ) 3 /Lo 




3.0 x 10 6 W kg" 


(CME) mass M cme 


m p n (Lo) 3 




5.7 x 10 13 kg 



A Parker-like solution for the solar wind is per construc- 
tion isothermal, i. e. 7 = 1, resulting in an adiabatic cooling 
for 7 > 1. In reality, the decrease in temperature T = p/p 
due to this adiabatic cooling of the expanding plasma is com- 
pensated by processes such as reconnective energy release 
and Alfvenic wave heating. A realistic inclusion of such ef- 
fects, while certainly desirable, is beyond the scope of this 
first approach, and, thus, reserved for future refinements of 
our model. As an alternative, we employ an ad hoc heating 
function 



Q = T C a(r) (V • u) 



(8) 



with a prescribed heating profile a(r) and a target tempera- 
ture T c . To derive the form of Eq. ([8]), we first seek the heat- 
ing function Qi SO which maintains a constant temperature T c 
everywhere, irrespectively of 7. This is done by inserting the 
corresponding isothermal equation of state 



P = T c p 



(9) 
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into the MHD equations ([T}j4]) and solving them analytically 
for Q, which yields 



Q|t=t c = (V >u)T c 



(10) 



Therefore, local heating (or cooling) at any heliocentric ra- 
dius r*h is conveniently achieved by choosing a(r*h) > 1 (or 
a(rh) < 1) in Eq. ffih. Test runs demonstrating the validity 



of this method have been carried out by |Kleimann| \ 2005 ). 



3 Numerical implementation 

3.1 Algorithm 

In order to integrate Eqs. ([Ifj4]) forward in time, we em- 
ploy a 3-D variant (Kleimann et al.| |2004j ) of a recent 
semi-discrete Central weighted essentially non-oscillatory 
(CWENO) scheme by |Kurganov and Levy| ( [2000| ) with third- 
order Runge-Kutta time stepping. Notable advantages of 
CWENO include its third-order accuracy in smooth regions 
(which automatically becomes second order near strong gra- 
dients to minimize spurious oscillations) and an easy general- 
ization to multi-dimensional systems of equations due to the 
fact that no (exact or approximate) Riemann solver is needed. 
The CWENO scheme thus allows simultaneously to achieve 
high shock resolution comparable with the best shock cap- 
turing schemes and high order convergence in smooth re- 
gions dominated by plasma waves. Although this scheme 
is not strictly total variation diminishing (TVD), simulations 



by |Levy et aL] ( |2000| ) do indicate an upper bound for the to- 
tal variation of their solutions. Moreover, |Havlfk and Liska| 
2006 ) use a set of astrophysically relevant test cases to com- 



pare the performance of several methods for ideal MHD and 
stress CWENO 's superior accuracy. 

3.2 Example test case: Alfven wings 

Various elementary tests of our implementation have been 
completed successfully ( [Kleimann et al.| 2004| ), such as ad- 
vection in one and two dimensions (also for propagation di- 
rections inclined at angles < v < tt/2 to the coordinate 
surfaces), shock tubes with and without magnetic field in- 
clusion, etc. 

While those standard tests will not be reproduced here, one 
rather advanced test setting, which is also of astrophysical 
relevance involving so-called "Alfven wings" is worth being 
mentioned. While the the finite extent of the wave-generating 
obstacle does not allow for an exact analytical solution, the 
usefulness of this simple but meaningful test problem stems 
from the fact that it incorporates several types of MHD waves 
(Alfven and slow/ fast magnetosonic), the expansion speed 
and characteristics of which can be verified quantitatively 
with theoretical expectations to ensure proper implementa- 
tion of the relevant physics. As shown by Drel l et aL] ( [1965| ), 
the movement of a conductive obstacle (e.g. a satellite or 



small planet) through a homogeneous fluid with a perpendic- 
ular magnetic field will generate standing MHD waves in the 
(u,B) plane called Alfven wings. This phenomenon plays a 
major role for the interaction of the moon Io with the Jovian 
magnetosphere (Neubauer, 1980; Linke r^eFaLl [1988 ), and 
also for artificial satellites in the magnetosphere of Earth, see 



e. g. Kopp and Schroer ( 1998) and references therein. 
The corresponding numerical test case, which works well 
both in 2-D and 3-D, involves an initially homogeneous flow 
u = uo e x of constant density, which is combined with a per- 
pendicular, equally homogeneous magnetic field B = B e z . 
A solid, spherical obstacle is then implemented by perform- 
ing an artificial deceleration 



u(r, t) u(r, t) x [1 — min(t, 1)] 

x [1 - tanh(4 max(||r|| - 1,0))] 



(11) 



after each time step, such that for t > 1, the flow will van- 
ish within ||r|| < 1. Figures [T] and [2] illustrate the emanat- 
ing wing structure. Direction and speed of propagation agree 
well with their respective theoretical expectations. 




Fig. 1. 3-D structure of a pair of Alfven wings, illustrated as an 
isocontour plot of absolute velocity. 



3.3 Choice of coordinates 

At first sight, the Sun's obviously spherical shape would sug- 
gest the use of spherical coordinates [r, tp], especially since 
the radial convergence of lines of constant (p entails the ad- 
ditional benefit of increased spatial resolution near the Sun's 
surface. On the other hand, the Courant-Friedrichs-Lewy 
(CFL) criterion of numerical validity and stability ( [Courant 
eTaLj [1928) , which requires the "velocity" Ax/ At to be 
greater than the maximum physical propagation velocity, im- 
poses a limit on the time step At based on the cell size Ax. 
The very choice of a coordinate system with varying grid cell 
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Fig. 2. Selected magnetic field lines in the (ix, B) plane, with color 
denoting flow velocity in normalized units. The flow is incident 
from the left. Since sound speed and Alfven speed are both unity, 
the wings emanate in a tailward 45° wedge when viewed in the 
obstacle's rest frame. The central circle marks the spherical volume 
inside of which deceleration according to Eq. flT) is applied. 



sizes, together with the requirement that the time step be uni- 
form on the entire grid, thus implies that At will be set by 
the Ax of the grid's smallest cell. For spherical coordinates, 
this means that the increased resolution at small r, however 
welcomed for physical reasons, would force At to be much 
lower than what the CFL criterion would require for most 
parts of the computational domain. This 'problem of small 
time steps' is avoided by using Cartesian coordinates, which 
have equal cell spacing everywhere and thus do not waste 
computing time on the larger cells. Even worse is the prob- 
lem of coordinate singularities at the poles $ G {0, 7r}, which 
require delicate numerical treatment. For these reasons, we 
opt for Cartesian coordinates [x,y, z], for which numerics 
are faster, simpler (esp. with respect to multi-dimensional ex- 
tension), and more stable. This is especially true since our 
CWENO code is built within a framework that allows for 



Cartesian Adaptive Mesh Refinement (AMR, see ( Kleimann 



et al. 2004|)). This is of high interest for more detailed stud- 



ies of, e.g., the inner structure of a CME. While AMR can, 
in principle, be used with spherical coordinates, its advan- 
tage is over-compensated by the fact that the convergence of 
grid spacing implies unacceptably low CFL numbers. (Note 
also that since CMEs generally do not exhibit any clear spa- 
tial symmetry, the use of non-Cartesian coordinates is not 
expected to entail any particular advantages for their descrip- 
tion.) 

3.4 Divergence cleaning 

Like many other algorithms, CWENO does not exactly 
conserve the solenoidality condition V • B = for the 



magnetic field, and a correction scheme becomes mandatory 
to avoid unphysical artifacts. From the wealth of existing 
schemes (for an overview see, e. g., Toth (2000)), we have 
evaluated the performance of the Generalized Lagrange 
Multiplier (GLM) approach by pedner et al.| ( |2002] ) against 
a classical projection scheme (see Sect. 3.4.2|). 



3.4.1 The GLM scheme 

The GLM scheme solves an additional equation 

d t ^ + (v f ) 2 V-B = -(v f /A) 



(12) 



for a position- and time-dependent Lagrange multiplier 
and adds a term — V\£ to the right hand side of Eq. ([3]). This 
procedure causes ^ (and hence V • B) to be damped with 
decay constant ra := A/^f, while at the same time advection 
of \I/ towards the boundary of the computational volume oc- 
curs at the highest permissible speed Vf (chosen to equal the 
global maximum of the fast magnetosonic speed in this case). 
Following |Dedner et al.| ( |2002] ), a value of 0. 1 8 is used for the 
second constant A. 

The main advantage of this method is that Eq. ( [T2| ) already 
possesses the correct conservative form, allowing for direct 
treatment with CWENO. In particular, physical conservation 
laws are not affected in any way. 

Figure [3] compares the performance of the two methods for a 
standard run. The obviously inferior performance of GLM 
can be explained by the fact that within a spherical layer 
C around the inner (solar) boundary, the boundary proce- 
dure described in section[5J]entails an averaging of the inner 
boundary value B in and the newly computed outer solution 
B out via 



^avg •— fBin + (1 



f)B a 



(13) 



for some function /:£>->■ [0, 1], which is bound to intro- 
duce a marked violation of the divergence constraint due to 
the first term of 



V/ • (Si, 



fV-B 



(14) 



being clearly non-zero. This divergence-laden field is then 
advected outwards by the wind flow, thus causing the mag- 
netic field to quickly become non-solenoidal in the outer re- 
gion as well. (This behavior becomes particularly evident in 
the left plot of Fig. [3]) 

Since the resulting magnitude of V • B in the non-solenoidal 
interface layer is inversely proportional to the layer's thick- 
ness, the problem cannot be avoided by choosing a differ- 
ent matching method (i. e. a different matching function 
r 4 /(r)). Note that this line of reasoning includes the case 
of doing no averaging at all: This simply corresponds to the 
limiting case / = / step , where 



step 



: r i — y 



1 : H 
0: IHI 



< 1 
> 1 



(15) 
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Since this non-solenoidal layer is actively re-created every 
time the newly computed outer solution is connected to the 
inner boundary, we may conclude that a suitable divergence 
cleaning procedure must kill the divergence immediately 
afterwards in one step (as the projection scheme does), 
rather than only damping/ transporting it away on somewhat 
longer timescales (GLM). 

We must therefore conclude that for investigations of this 
kind, the presence of an inner inflow boundary is, at least, 
difficult and may, in some cases, even preclude the use of 
the GLM scheme for divergence cleaning. (Note however, 
that the applicability of GLM to other settings lacking such 
an internal boundary remains unimpeded by this finding.) 



3.4.2 The projection scheme 

The so-called 'projection method' was originally developed 
by Chorin (1967]) for simulations of inviscid flow, and later 
applied in the context of MHD simulations by Brackbill and 
Barnes (1980). It solves the Poisson equation 



V 2 & = V'B 



(16) 



for and then subtracts V3> from B to ensure V • B = 0. 
While numerically expensive, it is able to reduce divergence 
errors down to machine accuracy, and will therefore be used 
in all simulations presented here. 



4 Boundary and initial conditions 

4. 1 Types of boundaries 

The computational volume consists of a brick-shaped region 
of space covering 100 x 70 x 50 cells in the x, y, and z di- 
rection, respectively. Each cell is a cube with a side length 
of typically Ax = Ay = Az = 0.1, implying a coverage of 
[5, 7, 10] Rq of real space. (We note that this relatively 
coarse resolution was chosen deliberately do demonstrate the 
excellent symmetry-maintaining properties of the employed 
scheme, see also Fig. [5] of Sect. 52_ Higher spatial resolu- 
tion, however desirable for the study of fine-scale structures, 
would tend to diminish the magnitude of numerical artifacts 
by which the scheme's performance could be judged, thereby 
hampering the usefulness of this demonstration.) The Sun's 
center is located at the origin, with the dipolar axis pointing 
into the positive z direction. The computational domain is 
surrounded by two layers of 'ghost cells', whose values are 
updated after each time step either from symmetry consid- 
erations (for 'mirror' boundaries intersecting the origin), or 
use of outflowing boundary conditions (at the actual 'outer' 
boundaries). 

The solar surface, which is represented by a sphere of unit 
radius located inside the computational volume, obviously 



does not coincide with any of the Cartesian coordinate sur- 
faces, and therefore requires special treatment, which is dis- 
cussed in detail in Sect. |5.1| and |5.2| This inner boundary 
is particularly delicate since it constitutes the surface from 
which the solar wind emanates, such that numerical artifacts 
imposed by an imperfect treatment of this boundary will be 
quickly advected through the entire domain. 

4.2 Initial conditions 

The generic setup for quiet-Sun solar wind simulations is as 
follows: At t = 0, the simulation is initialized with a radially 
symmetric wind flow u(r) = u(r) e r with 



u(r) 




(17) 



such that a super-sonic value u m is reached at the innermost 
boundary point, r m . This is done to ensure that the initial 
velocity at the outer boundary is as small as possible to al- 
low large time steps At, while at the same time being large 
enough to prevent numerical boundary artifacts from being 
transported inwards. 

The density scales as p(r) oc 1/r 3 , and the temperature is 
equal to a constant T c . The initial magnetic field is imple- 
mented using 



|t=o 



'Fir) 
V x [ — ^ sin i? e u 



(18) 



where F(r) = Po/r yields a dipole of strength Po that is 
aligned with the z axis. 

Since the projection scheme described in Sect. |3.4.2| w ill 
operate on the entire grid, the singularity of Eq. (|18|) at 
r = must be avoided. This is achieved by choosing a 
suitably matched polynomial for F(r) inside some small 
sphere around the origin. (Note that the radius of this sphere 
must be chosen at least several grid cell sizes smaller than 
unity to prevent the non-zero current density associated with 
F(r) ^ Po/r from causing unphysical Lorentz force accel- 
erations just outside the r = 1 boundary.) 

5 Numerical treatment of the solar surface boundary 

5 . 1 The interpolation method 

The inner (solar surface) boundary, which is just the sphere 
S := {r\ \\r\\ = 1}, obviously does not coincide with any of 
the Cartesian coordinate planes, which brings up the question 
of how these boundary conditions are best represented on the 
grid. Simple-minded attempts, such as keeping cell values 
inside the Sun fixed and integrating only those outside, have 
been tried but were shown to result in block-like artifacts 
at small radii (essentially tracing the envelope of the set of 
grid cells considered 'inside') where the problem's symme- 
try would stipulate spherical contours. While these artifacts 



J. Kleimann et al.: A novel MHD code for CME expansion 



div(B)/ Igrad ( I B I) I , time = 1.50 



div(B)/lgrad(lBl) , time = 1.50 





Fig. 3. Normalized divergence error n := (V • B)/\\ VVB • B\\ in the (poloidal) (x, z) plane for a standard solar wind run at time t = 1.5 
without correction (left) and using GLM (right). The improvement is substantial but still insufficient due to massive divergence values 
introduced at the inner boundary (unit circle around the origin). The projection scheme achieves k ~ 10 -6 (not shown). Note the different 
color scales. 



would of course diminish as spatial resolution is increased, 
it seems vital to obtain a high degree of symmetry-keeping 
already at this relatively coarse resolution, especially in view 
of the high numerical costs associated with increasing the 
number of grid cells in a 3-D simulation. 
After several possibilities have been tried, the following pro- 
cedure was adopted: 

1 . At initialization, all grid points which are located out- 
side S but have at least one of their 3 3 — 1 = 26 neigh- 
bors inside S are stored in a list X of 'interface points'. 
(The set neighbors of a cell is defined as the set of 
cells TV?/*,, with \i - z'|, \j - \k - k'\ G {0, 1} ex- 
cluding rijk itself.) 

2. After each time step (which only advances grid points 
outside S in time), a weighted average for each variable 

w G {p, pu x , pu y , pu z , B x , By, B z ,e} is computed for 
each n G X via 



with 



Y^iL^)- 1 w(r' a )\ I ($>,„) 



(19) 



{ 



77 r a H 



: r c 
S:r f 



outside of S 
inside of S 



(20) 



and Lj a := \\n — r' a ||, where the sums in Eq. ( 19 ) are 



taken over all neighbors of rj. The choice of weights 
ex (L/q,) -1 ensures that for ||r/|| — » 1, wj smoothly 
tends to the appropriate boundary value. Figure [4] serves 



to illustrate the situation. 

When the above procedure is applied to the Cartesian 
components of vectors such as u and B, it will usually 
destroy any possibly existing symmetry of these vector 
fields (e. g., if u is purely radial, the averaged u vec- 
tors will slightly deviate from the radial direction). In 
order to preserve such symmetries, all Cartesian vector 
components entering the averaging process of Eq. ( [T9| ) 
are first rotated until they are parallel to r/ before the 
averaging takes place, thus ensuring that the symmetry 
is preserved. 

3. In order to guarantee that the newly computed grid val- 
ues for X are independent of the ordering within that list, 
all computed averages are first stored in a separate field. 
Only when all the wj are known will they be copied 
onto the actual grid. 

Note that step 1 is executed only once, while steps 2 and 3 
are called after each integration time step. 
The above procedure gives the best results when applied to 
a scalar field that varies approximately linear in space. Near 
the solar surface, however, strong radial gradients of density 
are present. For this reason, it has been found to be advan- 
tageous to artificially reduce the density gradient in Eq. ( [T9] ) 
by multiplying p(r' a ) with ||r , ^|| 3 "" 4 before averaging, and 
consequently dividing pj by ||r/|| 3 -- 4 afterwards. 

5.2 Velocity extrapolation versus fixed boundary 



The averaging procedure of Sect. |5.1| keeps all quantities 
fixed on S. However, if solar wind configurations such as the 
Parker wind solution (Parker 1958| ) are to be reproduced, it 
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Fig. 4. 2-D analog of the averaging procedure. At each grid point 
ri Gl (shaded central box), a weighted average is computed from 
time-advanced values taken at neighbors outside S (crosses), and 
boundary values taken in the direction of neighbors inside S (cir- 
cles). The factors Li a entering into Eq. (l9\ are equivalent to the 
length of arrows in the diagram. 



seems questionable to apply this procedure to the velocity, 
since the requirement that r \-> must pass through 

a critical (sonic) point completely determines the solution 
topology, and thus eliminates the freedom to prescribe a fixed 
(Dirichlet) boundary value at r = 1. 

Different possibilities are conceivable to handle this prob- 
lem: 

1. Allow the velocity near S to adjust freely by radial in- 



ward extrapolation of the time-advanced solution ( Kep 



pens and Goedbloed||2000|), or 



2. enforce a fixed value for u on S in spite of the above 
problem, and accepting that (hopefully small) inaccura- 
cies will be introduced at small radii (Manche ster et al.[ 
[2004] ). 

While the second alternative is just what the above averaging 
procedure does, the first option, while being straightforward 
in spherical coordinates, is clearly non-trivial to implement 
on the present Cartesian grid. 

In analogy to the averaging scheme used for the other vari- 
ables, the adopted procedure (which replaces the procedure 
of Sect. |5.1| for u) is as follows: 

1. Prior to initialization, a list A of all grid points r a 
with 1 < \\r a || < 0.5 is set up and sorted by decreasing 
1 1 1*^ || (such that the outermost points will be processed 
first). 

2. For each va £ A, a sub-list of grid points r a,% is cre- 
ated, such that 



(i) ||r A || < \\vaA\ and 

(ii) \\r A - r A ,i\\ < r (with r « (2.. .3) Ax). 

In other words, the sub-list for r a contains grid points 
close to ta which are located at larger radii than r a 
itself. (Note again that steps 1 and 2 are executed only 
once.) 



3. After each time step, the radial mass flux 

f A ,i := (pu) A ,i ■ r A ,i \\rA,i\\ 



(21) 



is computed from the sub-list at each va, 
and a least- squares fit of the linear func- 
tion gA '• t i— >• Co, a + £\,a t is used to find 
the mass flux at ta (which is then given by 
gA(\\r A \\) = c ,a + Ci, A ||r A ||). Finally, the corre- 
sponding radial momentum is immediately afterwards 
written to the grid, such that its value is available to the 
extrapolation at the next point in the list. 

Figure [5] shows a comparison of both methods for the (un- 
magnetized) Parker wind case. The interpolation method's 
superior performance is in the range of a few percent only 
and has to be gauged against its increased computational ef- 
fort. Consequently, all of the simulations presented here use 
the Dirichlet method. (We note, however, that this may not 
always be appropriate when different parameter ranges are 
used. For instance, a higher base temperature T c will move 
the sonic point sunwards, leading to a higher flow velocity at 
the solar surface, and a presumably larger discrepancy to the 
zero- velocity condition.) 



6 Solar wind and CME simulations 

6.1 Creating equilibrium wind solutions 

For the initialization of our CME expansion studies, we first 
seek a well-defined MHD equilibrium resembling a 'quiet' 
(i.e. stationary) setting during solar minimum. While this is 
of course not strictly required for such studies, it is neverthe- 
less vital for the interpretation of the obtained results, since 
only then can structures like CMEs be clearly disentangled 
from the dynamics of the background flow. 
For magnetized, isothermal (7 = 1) winds, the system starts 



from the initial conditions of Sect. 4.2 and then quickly 
(within a few sound crossing times) settles into a stable equi- 
librium similar to the one depicted in the first frame of Fig. [6] 
At a distance of 5 Rq from the origin, the outflow velocity in 
x direction differs from that at the polar field line by a factor 
of about 

||u(5,0,0) 



1.26 c s 



160 km/s 



0.4 



||u(0,0,5)|| 3.13 c s 400 km/s 

which is due to the retaining force of the closed magnetic 
field lines in the equatorial region, and reminiscent of the 
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Fig. 5. Comparison of two possible methods to impose a boundary condition for u at the inner boundary r — 1: Dirichlet boundary 
condition (left) enforcing u\ r =\ — versus inward extrapolation (right). Both diagrams show scatter plots of absolute velocity versus 
heliocentric radius for all 40 3 grid points as the system converges towards the isothermal Parker wind solution (black solid curve). Colors 
are used to denote the moment of initialization (t = 0, blue), an intermediate step (t = 1, green), and the situation after the stationary 
equilibrium has been reached (t = 8, red). Since all grid points are shown, the scatter at a given radius can be seen as a measure of the 
simulation's spurious departure from radial symmetry. 



speed difference between the fast and slow solar wind. Ex- 
amples of non-isothermal hydrodynamical runs integrating 
the full energy equation ^ with the heating source term ([3} 
can be found in (Kleimann, 2005). 

It is noteworthy that essential features of the quiet inner he- 
liosphere, such as a latitudinal dependence of outflow veloc- 
ities resembling the fast and slow solar wind and the pole- 
ward transition from closed magnetic field lines (which span 
a static 'dead zone') below about 40° of latitude to an open, 
more radial field, are self-consistently reproduced by our 
model. In particular, it was found to be unnecessary to in- 
voke the method of latitude-dependent inner boundary condi- 
tions used by other authors ( Keppens and Goedbloed 2000 ; 



Man chester et al.| |2004| ) to reproduce this dichotomy: The 
magnetic dipole strength P proved fully sufficient to control 
the latitudinal extent of the closed-field helmet zone. As can 
be intuitively expected, a stronger B field at the surface will 
tend to conserve its arch-shaped closed structure, while in the 
limit Po —> 0, all field lines will be stretched out radially by 
the flow, and spherical symmetry is recovered. The choice of 
Po = 4 results in the intermediate case with an open/ closed 
transition near ±40° of solar latitude. 



6.2 Initialization of CME onset 

The present investigation focuses on the aspect of CME prop- 
agation, rather than on their actual nascency. Therefore, a 
simplifying approach similar to the one already employed 
by |Groth et~aL] ( |2000| ) and |Keppens and Goedbloed| (2000) 
will be used. This approach is based on a time-dependent 
boundary condition at the solar surface, generating a tran- 
sient, isothermal increase in density (and thus in pressure). If 
chosen sufficiently strong, this density excess is able to tear 
open the equatorial helmet streamer, causing the detachment 
of the excess matter as a rapidly expanding bubble. 
In order to initiate an eruption in the time interval 



°cmei ^cme 



T := [t c 

an additional, localized mass flux u a dd with 



Padd(V,£)|r= 
^addO,£)|r= 



Pcme V 
^cme &r 



(22) 



(23) 



is released at a pre-defined location on the solar surface (im- 
plying ||r|| = 1 for the remainder of this section). Without 
loss of generality, let the center of the eruption region be in 
the plane cp = 0, such that its location is just 



(24) 





•^cme ' 




f sin cme \ 


r - | 

' cme • — 1 


Verne 


l-l 


° 




\ ^cme / 




I cos$ cme / 
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For fixed time t, the value of p C me{r, t) should only depend 
on the angular distance 

ot(r, r cme ) = arccos(r • r cme ) (25) 
between r and r cme , such that p cme (r,t) possesses axial 



symmetry with respect to the r cme axis. Following Keppens 



and Goedbloed ( 2000 ), we employ the function 

fo E(r,t) :teT Aa<5 

c 

Pcme iy -> t) '•- 

with 

E(r,t) := sin 2 



: else 



t-t c 



cos 



2 / 7T a(r,r cme ) 



(26) 



(27) 



which connects smoothly to the undisturbed state in both 
space and time. Here 2S cme denotes the angular diameter of 
the circular eruption region ft cme , which is defined as the re- 
gion where p C me( r ? i) gives a non-zero contribution accord- 
ing to Eq. ([26]), and which thus covers a total solid angle 



(28) 



y cme 

J 2tt sin a da = 2tt [1 — cos(S c 



on the Sun's surface. The total mass released by the CME's 
eruption can be estimated as 



M CA 



3 (r,t) u cme duj dt 



(29) 



r a 
fo u, 



7T 2(S cme ) 2 - 7T 2 (1 - CQs£ cme ) 
CmeT " Cme 2 (^cme) 2 -^ 



with being the solid angle. For S cme = 30° = tt/6, this 
translates to physical units as 

MJO ^cme ^cme ^14 n /om 

cme,phys. ~ ^ X 10 i4 kg , (30) 

a typical value for a strong CME. 
6.3 CME expansion runs 

CME expansion runs have been carried out at various combi- 
nations of CME strength, heliographic latitude, dipolar field 
strength, etc. We first describe typical runs involving only a 
single CME, while the case of multiple events is deferred to 
the ensuing section. Unless indicated otherwise, the launch 
parameters fo = 16 and r cme = 1 were used. 



6.3.1 Single-event runs 

The panel of Fig. [6] shows selected snapshots of a typical 
simulation run involving an isolated CME event. The first 
frame depicts the equilibrium situation of the pre-emptive 
state. Following its initiation at the solar equator, the CME 



rapidly expands outwards, thereby quickly gaining both in 
size and speed. Note in particular the structures indicated 
by the kinked magnetic field lines that could lead to shocks 
in the non-isothermal case. In the third frame, while still 
continuing to accelerate, the CME reaches the volume 
boundary, thereby dragging the field lines outwards and 
deforming them almost radially. In the final frame, the CME 
has completely left the simulation volume and the system 
has relaxed into a state similar to the quiet initial situation. 

Taking advantage of the fully 3-D nature of our simula- 
tions, we can also access the dynamics in the perpendicular 
planes. Figure [7] shows time frames from the same run, 
this time viewed as a contour of the sharp, wall-shaped 
B z signature which arises when the CME runs into the 
background magnetic field and forces it to pile up ahead of 
it. As can be expected, the magnetic front moves fastest in 
the x direction, thus forming an elongated shell around the 
CME's core. On the opposite side, the CME's wake shows 
a marked reduction in field strength, which even includes 
an expanding region of reversed field direction trailing the 
CME. The region's growing extent is particularly evident 
from the dotted wedge discernible in Fig. [5] Note that the 
steep outward slope of ||£?|| (ex r -3 for a dipole) makes it 
necessary to normalize the values appropriately. 

To analyses the dynamics of the CME as a whole, a 
reliable tracer of its position is required. While the CME's 
density shows relatively large and irregular fluctuations 
which make it difficult to use it to monitor its location, 
we found the magnetic field signature of Fig. [7] to be more 
suitable for this purpose. Figure[8]may serve to illustrate this 
idea. From the resulting [t, x(t)] curves, we derive terminal 
velocities of 3.5 c s « 450 km/s and 4.5 c s « 580 km/s for 
CMEs launched with an initial velocity of u cme = and 2, 
respectively. 



6.3.2 Interacting CMEs 

With the rate of CME occurrence reaching several events 
per day during solar maximum, it is not unusual to find more 
than one CME to be present in a given section of interplan- 
etary space, a fact which motivates the numerical study of 
the interaction of CMEs. Simulations of this kind have been 



carried out by various authors (Vandas et al., 1997; Odstrcil 
[eTaL] |2003l |Schmidt and Cargill| |2004| |Wang et al.| |2005^ 7 



Interacting CMEs have also been linked to the modulation 
of type II radio bursts (Nunes, 2007 ), and their importance 
for SEP generation has been investigated by |Gop alswamy 
eTaL] p005] ) and |Vandas and Odstrcil| ( |200?|) u sing 2.5-D 



flux rope simulations. More recently, Lugaz| (|2008|) has 



connected earlier simulations ( |Lugaz et al. 2005, 



employing the BATS-R-US code ([Manchester et al. 



2007) 
2004) 



to actual LASCO data by means of synthetic observations. 
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Fig. 6. Time sequence of simulated CME expansion from pre-emption (t = 10.0) to expansion and return to equilibrium (at t = 20.0), 
showing contours of ||u|| (top) and log 10 n (bottom) in the poloidal plane (y = 0), with magnetic field lines superimposed. The CME speed 
at onset was chosen to be u C me — 2. An MPEG movie of this simulation, which covers the entire simulation from initialization (t = 0) to 
convergence into steady-state (near t — 10), CME expansion and back to near-equilibrium (t w 20), is available from the supplementary 
material at http://www.ann- geophys.net/0/l/2013/angeo-0- 1- 201 3- supplement. mpg 



While it is clear that at this initial stage, our simulations 
cannot be expected to rival the existing work in either detail 
or scientific content, we can nevertheless demonstrate our 
code's general applicability to this important sub-class of 
CME phenomenae. 



Figure [9] shows selected snapshots of a corresponding 
simulation run: At t = 10, a slow (u cme ,i = 0) CME 
is initiated along the x axis, to be quickly followed by a 
faster one (u C me,2 = 2) launched at t = 11.0 into the same 
direction. (Note that this terminology is merely used to 
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Fig. 7. Contour plot showing the value of B z in the (x, y) plane, normalized to the corresponding values at the pre-emptive equilibrium 
(t — 10.0). The unit circle marks the position of the solar surface. Note the field reversals occurring within the encircled regions (dotted), 
most notably in the CME's wake. 



distinguish both entities from each other. We do not intend 
to relate these to the slow/ fast dichotomy known from 
actual CME observations. As was shown at the end of the 
preceding section, both simulated CMEs would qualify as 
'slow' in this sense.) 

Both CMEs not only exhibit the individual effects of 
acceleration, expansion, and field line kinking already found 
and discussed in the previous case of an isolated CME, but 
there is apparently also a noticeable interaction between the 
two as the second CME gains speed and eventually collides 
with its predecessor. Note again the kinked magnetic field 
lines along with a corresponding density gradient, both 
related to discontinuities that would develop into shocks in 
the non-isothermal case. It is also interesting to observe 
that the prescribed density excess is sufficient to trigger a 
spontaneous, self-consistent outward acceleration, without 
the need to artificially 'push' the CME forward by enforcing 
a non-zero initial velocity at the instant of its launch. 



6.4 Connecting to observations 

Due to lack of in- situ data at small distances from the Sun, a 
direct comparison between simulation and actual CME data 
is currently not feasible. In order to at least qualitatively con- 
nect the simulations presented here to observations, five fixed 
locations at radii G {2, 4, 6, 8, 10} were chosen along the 
CME's trajectory (i. e. the x axis). At every time step, 
the values of the non- vanishing variables [B z , n,u x ] at these 
locations were extracted from the simulation data and then 



combined in the panel of Fig. 10 Thus, a time profile of 
these quantities is generated, as it would be seen by a sta- 
tionary observer while the CME moves past his location. (It 
should be noted that the profiles for particle density n and 
magnetic field B z at radius have been multiplied with 
(r^) 3 and — (r^) 2 , respectively, since otherwise the effect of 
radial dilution would not have allowed curves of various radii 
to be presented compactly in a single viewgraph. This obvi- 
ously only changes the relative size of two profiles against 
one another but leaves the shapes of individual profiles un- 
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Fig. 8. Height-time plot tracing the position of the normalized 
maximum of B z , which moves slightly ahead of the actual CME. 
Each vertical strip can be thought of as a cut along the x axis of 
Fig.[7]with identical color scale (including the dotted inversion line). 
The thick dashed line connects the respective maxima, thus forming 
an t H> x(t) position curve for the magnetic peak leading the CME. 
The solid line shows the corresponding x(t) plot for the faster CME 
(Ucme = 2 rather than 0). The corresponding contour stripes for this 
second CME are not shown. 



changed.) 

Using Fig. [TT] these plots can be contrasted with a compi- 
lation of the temporal evolution of the solar wind's MHD 
properties, as measured by the Advanced Composition Ex- 
plorer (ACE) for a magnetic cloud passing the probe's lo- 
cation, the inner Lagrange point at a heliocentric radius of 
0.99 AU. A number of qualitative similarities between ob- 
servation and our simulation can indeed be identified; espe- 
cially the sharp rise and slow decay of the velocity's maxima 
is clearly discernible in both cases. The sharp, almost needle- 
shaped peaks of the magnetic field profiles also exhibit a 
striking similarity. These pronounced field enhancements are 
induced by the CME's compression wave, and even seem to 
increase further as the driving CME accelerates outward. 
The notable differences in the total duration of passage 
(about one day for the magnetic cloud opposed to about one 
hour in the simulations) can easily be accounted for by the 
very different sites of observation. The cloud had much more 
time to extend from a presumably rather compact object to 
its full length of up to 1 AU. Also, the transit time cannot be 
expected to be totally independent of the duration of CME 
initiation (which in our case amounts to just 1.5 h real time). 
However, since observation and simulation stem from very 
different heliocentric radii, a direct, quantitative comparison 



between the respective profiles of Figs.[10|and[TT]is of course 
not feasible. Our attempts to identify common features be- 
tween them can therefore merely serve as a "reality check" 
on the general usefulness of these first simulation runs. Be- 



sides, they may serve to illustrate the type of comparison that 
are intended for future simulations covering the whole radial 
range up to Earth orbit. 



7 Conclusions 

We have reported on the creation of a 3-D MHD model of 
the near-Sun heliosphere, its numerical implementation and 
subsequent application to the propagation of CMEs. 
In order to adequately implement the Sun's spherical surface 
as an inner boundary on the Cartesian grid, a weighted av- 
eraging procedure was devised which is able to handle the 
huge gradients (most notably of mass density) present at this 
boundary. The use of this procedure also contributed to a re- 
duction of spurious departures from the problem's underly- 
ing symmetry, which result from the fact that the Sun's spher- 
ical (boundary) surface cannot be mapped to a Cartesian grid 
of finite cell spacing. Comparing a Dirichlet boundary con- 
dition for the velocity against free inward extrapolation, the 
latter was found to yield slightly more accurate results, albeit 
requiring a more complex numerical implementation. To en- 
sure a solenoidal magnetic field, the GLM scheme was found 
to be inappropriate due to the presence of internal bound- 
aries, and was thus abandoned in favor of a classical projec- 
tion method. 

After the model's CWENO-based numerical realization had 
satisfactorily passed various test cases, it was successfully 
employed to generate stable, self-consistent MHD equilibria 
of the quiet, magnetized solar wind. These were then them- 
selves used as initial configurations to simulate the expansion 
of CME-like plasma bubbles. Since the modelling is fully 
three-dimensional, the CME's direction of expansion can be 
chosen independently of the system's axis of symmetry; in 
particular, it is possible to study expansion within the eclip- 
tic plane. 

The extracted time profiles of density, flow velocity, and 
magnetic field strength show qualitative similarities to ac- 
tual in-situ data obtained from satellites at much larger helio- 
spheric distances. The fact that such similarities can be found 
lends support to the notion that the main physical processes 
which shape the structure of a CME occur shortly after on- 
set, whereas the ensuing phase of interplanetary propagation 
is merely characterized by dilution and (almost) self-similar 
expansion, although a direct simulation covering the entire 
range up to Earth orbit will be needed to make unambiguous 
statements about the CME's IP evolution and its persistent 
self- similarity (or lack thereof). In a future extension of this 
work, we intend to merge heated (i. e. non-isothermal) sce- 
narios with magnetized wind models, a step which, however 
desirable, could not yet be carried out due to remaining nu- 
merical difficulties. This direction seems even more promis- 
ing since both aspects have been proven to yield satisfactory 
solutions individually. 
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Fig. 9. Selected snapshots from a simulation of two interacting CMEs showing velocity (top) and decadic logarithm of density (bottom), as 
well as magnetic field lines (white) in the (y = 0) plane. The initial and final states are practically identical with the corresponding situation 
shown in the first and last frame of Fig. [6] and are therefore not repeated here. 



On the model side, we plan to include additional aspects into the CME's initialization to trigger its eruption. Since 
(such as localized heating and changes in magnetic topology) this will require a much higher grid resolution near the solar 
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Fig. 10. Simulated profiles of magnetic field, density, and fluid 
velocity as seen by static observers situated along the CME expan- 
sion direction at radii r / 'R® £ {2, 4, 6, 8, 10} for the sequence of 
Fig. [6] The minus sign at B z compensates for the magnetic field's 
north- south polarity (which has B z < at z = 0). Time is given in 
hours after CME onset. Temperature profiles are not shown due to 
7=1. 
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Fig. 11. Actual in situ measurements of a magnetic cloud near 1 
AU, adopted from Burlaga et al. ('2001). The general shape of these 
profiles is to be compared with the simulation results depicted in 
Fig.fTOl 



surface, a recourse to adaptive mesh refinement and/or paral- 
lelization becomes mandatory. 
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